APPROXIMATE VOLUME AND INTEGRATION FOR BASIC 
SEMI-ALGEBRAIC SETS* 
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Abstract. Given a basic compact semi-algebraic set K C M", we introduce a methodology 
that generates a sequence converging to the volume of K. This sequence is obtained from optimal 
values of a hierarchy of either semidefinite or linear programs. Not only the volume but also every 
finite vector of moments of the probability measure that is uniformly distributed on K can be 
approximated as closely as desired, and so permits to approximate the integral on K of any given 
polynomial; extension to integration against some weight functions is also provided. Finally, some 
Cn numerical issues associated with the algorithms involved are briefly discussed. 
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\^ ' 1. Introduction. Computing the volume and/or integrating on a subset 

K C M" is a challenging problem with many important applications. One possi- 
bility is to use basic Monte Carlo techniques that generate points uniformly in a 
box containing K and then count the proportion of points falling into K. To the 
best of our knowledge, all other approximate (deterministic or randomized) or exact 
techniques deal with polytopes or convex bodies only. Similarly, powerful cubature 
formulas exist for numerical integration against a weight function on simple sets (like 
^\J . e.g. simplex, box), but not for arbitrary semi-algebraic sets. 

^ ' The purpose of this paper is to introduce a deterministic technique that poten- 

7%^ i tially applies to any basic compact semi- algebraic set K C M". It is deterministic (no 

rp^ I randomization) and differs from previous ones in the literature essentially dedicated to 

fvj . convex bodies (and more particularly, convex polytopes). Indeed, one treats the orig- 

inal problem as an infinite dimensional optimization (and even linear programming 
(LP)) problem whose unknown is the Lebesgue measure on K. Next, by restricting to 
QQ ' finitely many of its moments, and using a certain characterization on the K-moment 

f^ . problem, one ends up in solving a hierarchy of semidefinite programming (SDP) prob- 

lems whose size is parametrized by the number of moments considered; the dual LP 
has a simple interpretation and from this viewpoint, convexity of K does not help 
/\ ' much. For a certain choice of the criterion to optimize, one obtains a monotone non 

^ . increasing sequence of upper bounds on the volume of K. Convergence to the exact 

value invokes results on the K-moment problem by Putinar [36. Importantly, there 
is no convexity and not even connectedness assumption on K, as this plays no role 
in the K-moment problem. Alternatively, using a different characterization of the 
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K-monient problem due to Krivine [23J, one may solve a hierarchy of LP (instead 
of SDP) problems whose size is also parametrized by the number of moments. Our 
contribution is a new addition to the already very long list of applications of the mo- 
ment approach (some of them described in e.g. Landau [24] and Lasserre [28]) and 
semidefinite programming [45]. In principle, the method also permits to approximate 
any finite number of moments of the uniform distribution on K, and so provides a 
means to approximate the integral of a polynomial on K. Extension to integration 
against a weight function is also proposed. 

Background. Computing or even approximating the volume of a convex body is 
hard theoretically and in practice as well. Even if K C M" is a convex polytope, 
exact computation of its volume or integration over K is a computational challenge. 
Computational complexity of these problems is discussed in e.g. BoUobas 7 and 
Dyer and Frieze [TT] . Any deterministic algorithm with polynomial time complexity 
that would compute an upper bound vol (K) and a lower bound vol (K) on vol (K) 

cannot yield an upper bound on the ratio vol (K)/vol(K) better than polynomial in 
the dimension n. Methods for exact volume computation use either triangulations 
or simplicial decompositions depending on whether the polytope has a half-space 
description or a vertex description. See e.g. Cohen and Hickey, [TO], Lasserre [25] . 
Lawrence [32] and see Biieler et al. p] for a comparison. Another set of methods 
which use generating functions are described in e.g. Barvinok ^ and Lasserre and 
Zcron [301. Concerning integration on simple sets (e.g. simplex, box) via cubature 
formulas, the interested reader is referred to Gautschi [T4l[T5] and Trefethren [43] . 

A convex body K C M" is a compact convex subset with nonempty interior. A 
strong separation oracle answers either x G K or a; ^ K, and in the latter case 
produces a hyperplane separating x from K. A negative result states that for every 
polynomial-time algorithm for computing the volume of a convex body K C M" 
given by a well-guaranteed separation oracle, there is a constant c > such that 
vol(K)/vol (K) < (en/ log n)" cannot be guaranteed for n > 2. However, Lovasz [33] 

proved that there is a polynomial-time algorithm that produces vol (K) and vol (K) 

satisfying vol (K)/vol (K) < n" (n + 1)"/^, whereas Elekes 13 proved that for < e < 

2 there is no polynomial-time algorithm that produces vol (K) and vol (K) satisfying 

vol (K) / vol (K) < (2 - e)". 

If one accepts randomized algorithms that fail with small probability, then the sit- 
uation is more favorable. Indeed, the celebrated Dyer, Frieze and Kanan probabilistic 
approximation algorithm [12] computes the volume to fixed arbitrary relative preci- 
sion e, in time polynomial in e~^. The latter algorithm uses approximation schemes 
based on rapidly mixing Markov chains and isoperimetric inequalities. See also hit- 
and-run algorithms for sampling points according to a given distribution, described 
in e.g. Belisle [S], Belisle et al. [B], and Smith [IT] 

Contribution. This paper is concerned with computing (or rather approximating) 
the volume of a compact basic semi-algebraic set K C M" defined by 

K := {.tgM" : gj(a;) > 0, j = l,...,m} (1.1) 

for some polynomials (gj)7Li C M[x]. Hence K is possibly non-convex and non- 
connected. Therefore, in view of the above discussion, this is quite a challenging 
problem. 
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(a) We present a numerical scheme that depends on a parameter p, a polynomial 
that is nonnegative on K (e.g. p = 1). For each parameter p, it provides converging 
approximations of moments of the measure uniformly supported on K (with mass 
equal to vol(K)). For the choice p = 1 one obtains a monotone non- increasing 
sequence of upper bounds that converges to vol(K). 

(b) The way we see the problem dates back to the 19th century pioneer work 
in the one-dimensional case by Chebyshev [9j, Markov [Ml and Stieltjes [^, where 
given n moments Sk — J t^f{t)dt, k — 0, ...,n — 1, and a < c < d < b, one 

wishes to approximate the integral J f{t)dt and analyzes asymptotics as n — > oo; 
characterizing feasible sequence (s^) is referred to as the Markov moment problem 
(and L-moment problem if in addition one requires < f < L for some scalar L). For 
an historical account on this problem as well as other developments, the interested 
reader is referred to e.g. Krcin [5J, Krein and Nuldelman '^, Karlin and Studden 
PP] and Putinar [37]. 

Our method combines a simple idea, easy to describe, with relatively recent pow- 
erful results on the K-momcnt problem described in e.g. [351 IMl SD]- It only re- 
quires knowledge of a set B (containing K) simple enough so that the moments of the 
Lebesgue measure on B can be obtained easily. For instance B := {x G R" : ||.t||p < a} 
with p = 2 (the scaled n-dimensional ball) or p = cxd (the scaled n-dimensional box) 
and a G R a given constant. Then computing vol(K) is equivalent to computing the 
mass of the Borel measure /i which is the restriction to K of the Lebesgue measure 
on B. This in turn is translated into an infinite dimensional LP problem P with 
parameter p (some polynomial nonnegative on K) and with the Borel measure fi as 
unknown. Then, from certain results on the K-moment problem and its dual theory 
of the representation of polynomials positive on K, problem P can be approximated 
by an appropriate hierarchy of semidefinite programs (SDP) whose size depends on 
the number d of moments of /i considered. One obtains approximations of the mo- 
ments of /i which converge to the exact value as d ^ oo. For the choice p = 1 of the 
parameter p, one even obtains an non-increasing sequence of upper bounds converging 
to vol(K). Asymptotic convergence is ensured by invoking results of Putinar [36j on 
the K-moment problem. Alternatively, one may replace the SDP hierarchy with an 
LP hierarchy and now invoke results of Krivine [23] for convergence. 

Interestingly, the dual of each SDP relaxation defines a strenghtening of P*, the 
LP dual of P, and highlights why the problem of computing the volume is difficult. 
Indeed, one has to approximate from above the function / (= p on K and on B \ K) 
by a polynomial h of bounded degree, so as to minimize the integral /^(ft. — f)dx. 
From this viewpoint, convexity of K plays no particular role and so, does not help 
much. 

(c) Let d e N be fixed, arbitrary. One obtains an approximation of the moments 
of degree up to d of the measure /i on K, as closely as desired. Therefore, this tech- 
nique also provides a sequence of approximations that converges to /„ qdx for any 
polynomial q of degree at most d (in contrast, Monte Carlo simulation is for a given 
q) . Finally, we also propose a similar approximation scheme for integrating a polyno- 
mial on K against a nonnegative weight function w{x). The only required data are 
moments of the measure dv — wdx on a simple set B (e.g. box or simplex) containing 
K, which can be obtained by usual cubature formulas for integration. 

On the practical side, at each step d of the hierarchy, the computational workload 
is that of solving an SDP problem of increasing size. In principle, this can be done 
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in time polynomial in the input size of the SDP problem, at given relative accuracy. 
However, in view of the present status and limitations of current SDP solvers, so far 
the method is restricted to problems of small dimension n if one wishes to obtain 
good approximations. The alternative LP hierarchy might be preferable for larger 
size problems, even if proved to be less efficient when used in other contexts where 
the moment approach applies, see e.g. [271131]. 

Preliminary results on simple problems for which vol (K) is known show that 
indeed convexity plays no particular role. In addition, as for interpolation problems, 
the choice of the basis of polynomials is crucial from the viewpoint of numerical 
precision. This is illustrated on a trival example on the real line where, as expected, 
the basis of Chebyshev polynomials is far better than the usual monomial basis. 
In fact, it is conjectured that trigonometric polynomials would be probably the best 
choice. Finally, the choice of the parameter p is also very important and unfortunately, 
the choice of p = 1 which guarantees a monotone convergence to vol (K) is not the 
best choice at all. Best results are obtained when p is negative outside K. 

So far, for convex polytopes, this method is certainly not competitive with exact 
specific methods as those described in e.g. i8_. It rather should be viewed as a 
relatively simple deterministic methodology that applies to a very general context for 
which even getting good bounds on vol (K) is very difficult, and for which the only 
alternative presently available seems to be brute force Monte Carlo. 

2. Notation, definitions and preliminary results. Let M.[x] be the ring of 

real polynomials in the variables x = (a;i, . . . , a;„), and let S^[a;] C R[a;] be the subset 
of sums of squares (SOS) polynomials. Denote M[a:]d C M.[x] be the set of polynomials 
of degree at most d, which forms a vector space of dimension s{d) — ("^ ). If 
/ G IR[x]d, write f{x) — J2aeN" fa^" in the usual canonical basis (x"), and denote 
by f = (fa) e W'''^^ its vector of coefficients. Similarly, denote by T,'^[x]d C S^[a:] the 
subset of SOS polynomials of degree at most 2d. 

Moment matrix. Let y = (i/a) be a sequence indexed in the canonical basis (a;") 
of M.[x], let Ly : M[x] ^ R be the linear functional 

/ (=^/„.t") ^ Ly(/) =5^/.2/o, 

a a 

and let Md{y) be the symmetric matrix with rows and columns indexed in the canon- 
ical basis (a;"), and defined by: 

Md{y){a,(i) := iy(x"+^) = y„+/3, 

for every a, /? £ N^ := {a G N" : |a| (= Y.^ "0 < d]. 

A sequence y — (j/q) is said to have a representing finite Borel measure ^ if 
Va = J x"dfi for every a G N". A necessary (but not sufficient) condition is that 
Md{y) h for every d G N. However, if in addition, \ya\ < M for some M and for 
every a G N", then y has a representing measure on [—1,1]". 

Localizing matrix. Similarly, with y = (y^) and g G M\x\ written as 



X ^ g{x) = Y^ 

7eN" 



9 J ■^ J 



let Md{gy) be the symmetric matrix with rows and columns indexed in the canonical 



Approximate volume and integration 5 

basis (x"), and defined by: 

Md{gy)ia,P) := Ly {g{x) x°'+'^) = ^g^y^+p+j, 

7 

for every a,/3 G NJJ- A necessary (but not sufficient) condition for y to have a 
representing measure with support contained in the level set {x : g{x) > 0} is that 
Md{gy) h for every d e N. 

2.1. Moment conditions and representation theorems. The following re- 
sults from the K-moment problem and its dual theory of polynomials positive on K 
provide the rationale behind the hierarchy of SDP relaxations introduced in [26j , and 
potential applications in many different contexts. See e.g. |28j and the many refer- 
ences therein. 

SOS-based representations. Let Q{g) C M[x] be the quadratic module generated 
by polynomials {gj)JLi C M.[x], that is, 



Q(g) ■■= \'^o + J2 "^ 9j ■■ {^j)T=i c s^N \ . (2.1) 

Assumption 2.1. The set K C M" m lll.l\) is compact and the quadratic poly- 
nomial X 1-^ a'^ — llxp belongs to Q{g) for some given constant a G R. 

Theorem 2.2 (Putinar's Positivstellensatz [36]). Let Assumption [2J\ hold. 

(a) // / e M[a;] is strictly positive on K, then f G Q{g)- That is: 

m 

/ = ao + ^a,g„ (2.2) 

for some SOS polynomials (crj)jLi C S^[a;]. 

(b) If y — iya) is such that for every d G N, 

Md{y)hO; Mdig,y)hO, j = l,...,m, (2.3) 

then y has a representing finite Borel measure fi supported on K. 

Given / G M[x], or y = (ya) C M, checking whether (|^ holds for SOS 



^2 

d fixed, reduces to solving an SDP, 



{(Jj) C ^"^[x] with a priori bounded degree, or checking whether (|2.3p holds with 



Another type of representation. Let K C B be as in (|l.ip and assume for 
simplicity that the gjS have been scaled to satisfy < gj < 1 on K, for every j = 
1, . . . , m. In addition, assume that the family of polynomials (1, 51, ... , gm) generates 
the algebra R[x]. For every a G N™, let 5" and (1 — g)^ denote the polynomials 

x^ gixr := giixr^ ■ ■ ■ gmixr^ 
and 

x^ (1 - g{x)f :- (1 - g,{x)f^ ■■■{!- g^{x)f-. 
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the following result is due to Krivine 23J but is explicit in e.g. Vasilescu |46j . 
Theorem 2.3. 

(a) If f G M.[x] is strictly positive on K, then 

/= E c-;3.g"(i-g)^ (2.4) 

q,/3gN™ 

for finitely many nonnegative scalars (cap) C IR+. 

(b) If y = {Ua) is such that 

Lyig'^ {1 - gf) > 0, (2.5) 

for every a,(3 d N™, then y has a representing finite Borel measure fi supported on 
K. 

Theorem 12.31 extends the well-known HausdorfF moment conditions on the hyper 
cube [0, 1]", as well as Handelman representation [17 for convex polytopes K C M". 
Observe that checking whether (|2.4p . resp. (|2.5p . holds with a, (3 bounded a priori, 
reduces to solving an LP in the variables (cq/j), resp. (ya)- 

2.2. A preliminary result. Given any two measures /ii, /i2 on a Borel cr-algebra 
S, the notation /ii < fi2 means /J.i(C) < /i2(C) for every C e S. 

Lemma 2.4. Let Assumvtion \2.1\ hold and let yi ~ (yi^) and yi = {y2a) ^c l-wo 
moment sequences with respective representing measures fii and /i2 on K. // 

Md{y2 - yi) ^ ; Md(.g, (ya - yi)) ^0, j = 1, . . . , m, 

for every d G N, then /ii < fi2- 

Proof As Md{y2 - yi) >= and M^gj (y2 - yi)) ^ for j = 1, . . . , m and d G N, 
by Theorem 12.21 the sequence yo := y2 — yi has a representing Borel measure /io on 
K. From j/Qq + 2/1q = 2/2q for every a £ N", we conclude that 



a;" d^o + / x" d^i = / x" d^2, Va G N" 

and as K is compact, by the Stone- Weierstrass theorem, 

/ d^o + / / d^i = / / d^J.2 



for every continuous function / on K, which in turn implies fio + fJ-i —1^2, i.e., the 
desired result /ii < /i2. □ 

3. Main result. We first introduce an infinite-dimensional LP problem P whose 
unique optimal solution is the restriction /i of the normalized Lebesgue measure on 
B (hence with /i(K) = vol (K)/2") and whose dual has a clear interpretation. We 
then define a hierarchy of SDP problems (alternatively, a hierarchy of LP problems) 
to approximate any finite sequence of moments of /i, as closely as desired. 

3.1. An infinite-dimensional linear program P. After possibly some nor- 
malization of the defining polynomials, assume with no loss of generality that K C 
B C [—1,1]" with B a set over which integration w.r.t. the Lebesgue measure is easy. 
For instance, B is the box [—1, 1]" or B is the euclidean unit ball. 
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Let B be the Borel cr-algebra of Borel subsets of B, and let /i2 be the Lebesgue 
measure on B, normalized so that 2"/X2(B) = vol(B). Therefore, if vol (C) denotes 
the n-dimensional volume oi C <E B, then H2(C) ~ vol (C)/2" for every C G B. 

Also, the notation fii <C /i2 means that /ii is absolutely continuous w.r.t. fX2, and 
Li(/i2) is the set of all functions integrable w.r.t. /i2. By the Radon-Nikodym theorem, 
there exists a nonnegative measurable function / G Li{^2) such that /Ui(C) = J^, /(^M2 
for every C S ,B, and / is called the Radon-Nikodym derivative of /ii w.r.t. 1x2- In 
particular, /ii < 112 obviously implies /ii <C /i2- For K G Z?, let M(K) be the set of 
finite Borel measures on K. 

Theorem 3.1. Let K G S with K C B and let p G R[x] be positive almost 
everywhere on K. Consider the following infinite- dimensional LP problem: 



P: sup { p d^i : ^ii < H2; A^i G M(K) } 
Ml J 



(3.1) 



with optimal value denoted supP (and maxP if the supremum is achieved). 

Then the restriction /i^ of fi2 to K is the unique optimal solution ofV and maxP = 
J pd^l ~ Jj^pdiJ,2- In particular, if p = 1 then maxP — vol(K)/2". 

Proof Let ^J be the restriction of ^2 to K (i.e. /zJ(C) = ^^2(0 n K), VC G S). 
Observe that nl is a feasible solution of P. Next, let ni be any feasible solution of P. 
As /xi < IJ.2 then 

Aii(cnK) <^2(cnK) = Mi(cnK), vcgS, 

and so, /ii < fxl because ni and nl are supported on K. Therefore, as p > on K, 
J pdfii < J pdfil which proves that /ij is an optimal solution of P. 

Next suppose that /xi ^ nl is another optimal solution of P. As ^1 < /i| then 
Hi <^ fj^l and so, by the Radon-Nikodym theorem, there exists a nonnegative measur- 
able function / G Li{fil) such that 

Mi(C) = f d^ii = f f{x)d^il{x), VCgSHK. 
Jc Jc 

Next, as jii < ji\, pL\ — jjLi =: fiQ is a finite Borel measure on K which satisfies 



< Mo(C) - ( {I- f{x))dpi\{x), VCGi3nK, 
Jc 

and so 1 > f{x) for almost all a; G K. But then, since J pd/ii — J pd/il, 



= / pdfio = / p{x){l- f{x))dHl{x), 

which (recalling p > almost everywhere on K) implies that f{x) — 1 for almost-all 
a; G K. And so fii — fi^. D 

3.2. The dual of P. Let T be the Banach space of continuous functions on 
B (equipped with the sup norm) and J-^ its positive cone, i.e., the elements f E J-^ 
which arc nonnegative on B. The dual of P reads: 

P* : M^ if fdti2 ■■ / > P on K} (3.2) 

with optimal value denoted inf P* (minP* is the infimum is achieved). 
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Hence, a minimizing sequence of P* aims at approximating from above the func- 
tion / (= p on K and on B \ K) by a sequence (fi) of continuous functions so as to 
minimize J fedfi2. 

Let X 1-^ d{x,K.) be the euchdean distance to the set K and with ei > 0, let 
Kf := {a; G B : d(x,K) < e^} be an open bounded outer approximation of K, so 
that B \ K^ is closed (hence compact) with e^ — > as £ ^ cx). By Urysohn's Lemma 
[U A4.2, p. 379], there exists a sequence {fi) C J-+ such that < /^ < 1 on B, /^ = 
on B \ Kf , and /^ = 1 on K. Therefore, 



f hdfi2 = yol{K)/T+ f fedfi2, 



KAK 

and so / fe.d^2 —* vol (K)/2" as ^ — * oo. Hence, for the choice of the parameter p = 1, 
vol(K)/2" is the optimal value of both P and P*. 

3.3. A hierarchy of semidefinite relzLxations for computing the volume 

of K. Let y2 = {y2a) be the sequence of all moments of fi2- For example, if B = 

[-1, 1]", then 

Let K be a compact semi-algebraic set as in ()l.ip and let rj — [(deggj)/2], j = 
1, . . . , TO. Let p 6 M.[x] be a given polynomial positive almost everywhere on K, and 
let To := [(degp)/2]. For d > maxj rj, consider the following semidefinite program: 

sup Lyi (p) 

s.t. Md(yi) h (3.3) 

Md(y2-yi) ^0 

Md-rj (5j yi) ^0, j = 1, . . . , TO 

with optimal value denoted sup Q^ (and max Q^ if the supremum is achieved) . 

Observe that supQd > maxP for every d. Indeed, the sequence y^ of moments 
of the Borel measure fil (restriction of fi2 to K and unique optimal solution of P) is 
a feasible solution of Q^ for every d. 

Theorem 3.2. Let Assumvtion \2. 1\ hold and consider the hierarchy of semidefi- 
nite programs (Qd) in \3.k^) . Then: 

(a) Qc; has an optimal solution (i.e. supQd = maxQ^j and 

maxQrf 4 / pdii2, as d —> oo. 

(b) Let yi'' ~ ivia) ^6 '^n optimal solution of Qd, then 

hm 2/1^ = / a;"d/^2, Va G N". (3.4) 

Proof, (a) and (b). Recall that B C [—1,1]". By definition oi fi2, observe that 
\y2a\ < 1 for every a G Nj,;, and from Md{y2 — Yi) h 0, the diagonal elements 
2/22Q ~ yi2a ^I'e nonnegative. Hence yi2a ^ ^220 for every a G N^ and therefore, 

max[yiQ, max Ly^{xf )] < 1. 

i—l. n 
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By [551 Lemma 1], this in turn implies that |yiQ,| < 1 for every a e ^2d' ^^"^ ^^ ^^^ 
feasible set of Q^ is closed, bounded, hence compact, which in turn implies that Q^^ 
is solvable (i.e., has an optimal solution). 

Let yi"^ be an optimal solution of Q^ and by completing with zeros, make yi'' an 
element of the unit ball Boo of loo (the Banach space of bounded sequences, equipped 
with the sup-norm). As Zoo is the topological dual of h, by the Banach- Alaoglu 
Theorem, Boo is weak • compact, and even weak * sequentially compact; see e.g. 
Ash [T]. Therefore, there exists yi* 6 Boo and a subsequence {dk} C N such that 
yj^'ifc -^ y^* as fc ^ oo, for the weak -k topology <7{loo, h)- In particular, 

lim yit = yi^, VaeN". (3.5) 

k — >oo 

Next let d e N be fixed, arbitrary. From the pointwise convergence p.Sp we also 
obtain Md{yi*) h and Md{y2 - Yi*) h 0. Similarly, Md^r^ig^yi*) h for every 
J = 1, . . . , m. As d was arbitrary, by Theorem 12. 2i yi* has a representing measure fj,i 
supported on K C B. In particular, from p.Sp . as fc — + oo, 

maxP < maxQdfc = L^d^{p) | Ly'(p) = / pdfii. 

Next, as both ^i and ^2 are supported on [—1, 1]", and Md{y2 — yi*) h for every 
d, one has \y2a — Vial < 1 for every a G N". Hence y2 — y^ has a representing 
measure on [—1,1]". As in the proof of Lemma l2.4n . we conclude that /^i < fj,2. 
Therefore fii is admissible for problem P, with value Ly«(p) = J pd^i > max P. 
Therefore, fii must be an optimal solution of P (hence unique) and by Theorem 13. 1[ 
Ly* = J pdfXi = J-^pdiJ,2- As the converging subsequence {dk} was arbitrary, it 
follows that in fact the whole sequence yi^ converges to yi* for the weak • topology 
<^{loo,li)- And so (|3.4p holds. This proves (a) and (b). D 

Writing Md(yi) = Ea^^yia: and Md-r^ig^yi) = J2a ^ivia for appropriate 
real symmetric matrices {Aa, B^), the dual of Qd reads: 

' ^inf^^ (Afd(y2),y) 

Qd- { s.t. (A„,r-x)-^(B^,Zj)=p„ 

X, Y, Zj y 0, 

where {X, Y) — trace (XY) is the standard inner product of real symmetric matrices, 
and X y stands for X is positive semidefinite. This can be reformulated as: 



Q^ 



inf / h (i/i2 

m 

s.t. h-p = ao + ^a.jgj (^'^^ 

h e s^Hd, (To e T.^[x]d, (Tj £ j:^[x]d^r,- 

The constraint of this semidefinite program states that the polynomial h—pis written 
in Putinar's form (|2.2p and so ft, — p > on K. In addition, h > because it is a sum 
of squares. 



^If K C [-1, 1]" then in Lemma [Ol the condition Ma{y2 - yi) h 0, Vd G N, is sufficient. 
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This interpretation of Q|J also shows why computing vol(K) is difficult. Indeed, 
when p = 1, to get a good upper bound on vol(K), one needs to obtain a good 
polynomial approximation h G R[x] of the indicator function I-k{x) on B. In general, 
high degree of h will be necessary to attenuate side effects on the boundary of B and 
K, a well-known issue in interpolation with polynomials. 

Proposition 3.3. //K and B\K have a nonempty interior, there is no duality 
gap, that is, both optimal values of Q^ and Q^ are equal. In addition, Q^ has an 
optimal solution {h* , (cr*)). 

Proof. Let fii be the uniform distribution on K, i.e., the restriction of /i2 to K, 
and let yi — (yi^) be its sequence of moments up to degree 2d. As K has nonempty 
interior, then clearly Md{yi) >- and Md-r {gj yi) ^ for every j = 1, . . . ,m. If 
B \ K also has nonempty interior then Md{y2 — yi) >- because with / G M.[x]d with 
coefficient vector f , 



(f,M,(y2-yi)f) = [ 

JB\K 



f{xfdii2. v/ e m.[x\d. 

\i 



Therefore Slater's condition holds for Qd and the result follows from a standard result 
of duality in semidefinite programming; see e.g. [iS]. D 

Remark 3.4. Let f G R[a;] and suppose that one wants to approximate the 
integral J* :— Jj^fd^2- Then for d sufficiently large, an optimal solution of Qd 
allows to approximate J* . Indeed, 

J* = fdn2 == fd^ii = Lyi-(/) = YI f^y^a^ 

where yi* is the moment sequence of fix, the unique optimal solution of P (the re- 
striction of fi2 to K.). And so, from |j'.4[ j, L d{f) « J* when d is sufficiently large. 

3.4. A hierarchy of linear programs. Let K C B C [— 1, 1]" be as in (|l.ip 
and assume for simplicity that the gjS have been scaled to satisfy < g^ < 1 on K for 
every j = 1, . . . ,m. In addition, assume that the family of polynomials (1, gi, . . . , (?„i) 
generates the algebra R[x]. For d G N, consider the following linear program: 



supy^ 2/10 

s-t. Ly,-yi ( 11(1 + x,r^ (1 - a;,)^* J > 0, a, /? G N^ 

LyAg'' {I - gf) > 0, a,/3GN2 



(3.7) 



with optimal value denoted supL^; (and maxhd if supL^; is finite). Notice that 
supLrf > vol(K)/2" for all d. Indeed, the sequence y| of moments of the Borel 
measure /i| (restriction of fi2 to K and unique optimal solution of P) is a feasible 
solution of hd for every d. 

Theorem 3.5. For the hierarchy of linear programs (L^) in |5'.?] ), the following 
holds: 

(a) Ld has an optimal solution (i.e. supL^ = maxL^J and maxL^ | vol(K)/2" 
as d —> oo. 

(b) Let yi'' be an optimal solution ofhd. Then iS.4\) holds. 
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Proof. We first prove that L^j has finite value. L^ always has a feasible solution 
yi, namely the moment vector associated with the Borel measure /ii, the restriction 
of fj,2 to K, and so supL^ > vol(K)/2". Next, from the constraint Ly2_yj(») > 
with a = (3 = 0, we obtain j/iq < 2/20 ^ 1- Hence supL^ < 1 and therefore, the 
linear program L^ has an optimal solution yi**. Fix 7 G N" and e > 0, arbitrary. As 
|a;T| <l<l + eonB (hence on K), by Theorem I^I^a). 

a,/3eN™ 

for some (c^a) C M+ with |a|, \(3\ < s-y. Hence, as soon as d > Sj, applying Ly^d 
yields 

(1 + e) yig ± yi^ = ^ c^^ Ly, (5" (1 - gf) > 0, 



and so 



V7eN": |2/4l<(l + e)2/i^<l + e, Vd > s^. (3.8) 



Complete yi"* with zeros to make it an element of E°°. Because of p.Sp . using a 
standard diagonal element, there exists a subsequence (d^) and an element yi* € 
(1 + e) Boo (where Boo is the unit ball of ^00) such that (|3.5p holds. Now with 
a,/3 e N™ fixed, arbitrary, dSH) yields Ly,,{g°' (1 - g)'^) > 0. Hence by Theorem 
I2.3r b). yi* has a representing measure m supported on K. Next, let yo '■— y2 — yi*- 
Again, (|3.5p yields: 

and so by Theorem I2.3r b) , yo is the moment vector of some Borel measure /iq sup- 
ported on [—1, 1]". As measures on compact sets are identified with their moments, 
and 2/Oq + Via = V^a fo'" every a S N", it follows that /xq + /ii = ^2, and so /ii < ^2- 
Therefore, /ii is an admissible solution to P with parameter p = 1, and with value 
/ii(K) = 2/ig > vol (K)/2". Hence, /ii is the unique optimal solution to P with value 
/^i(K) =vol(K)/2". 

Finally, by using p.Sp and following the same argument as in the proof of Theorem 
13.21 one obtains the desired result p.4p . D 

Remark 13.41 also applies to the LP relaxations (|3.7p . 

3.5. Integration against a weight function. With K C B as in (|l.ip suppose 
now that one wishes to approximate the integral 

J* ;= / f{x)w{x)dx, (3.9) 

for some given nonnegative weight function w : M" —^ M, and where / G ^[xjt; is some 
nonnegative polynomial. One makes the following assumption: 

Assumption 3.6. One knows the moments y2 = (2/2 q) of the Borel measure 
dfi2 = wdx on B, that is: 

y2a = / x°'dfi2 (^ f x"wix)dxj , aeN". (3.10) 
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Indeed, for many weight functions w, and given d e N, one may compute the moments 
y2 = (y2a) of fi2 via cubature formula, exact up to degree d. In practice, one only 
knows finitely many moments of fi2 , say up to degree d, fixed. 
Consider the hierarchy of semidefinite programs 



sup 


^yiU) 




yi 






s.t. 


Af<i(yi) 


h 




Mdiy2-yi) 


h 




Md-r, {g J yi) 


h 



(3.11) 



J-1, 



with y2 as in Assumption 13.61 

Theorem 3.7. Let Assumption \2. 1\ and \3.6\ hold and consider the hierarchy of 
semidefinite programs (Qd) in 113.11]) with y^ as in li3.1U\) . Then Q^ is solvable and 
maxQd I J* as d ^ CO. 

The proof is almost a verbatim copy of that of Theorem [ 



4. Numerical experiments and discussion. In this section we report some 
numerical experiments carried out with Matlab and the package GloptiPoly 3 for ma- 
nipulating and solving generalized problems of moments [18j . The SDP problems were 
solved with ScDuMi 1.1R3 [35]. Univariate Chebyshev polynomials were manipulated 
with the chebfun package [3]. 

The single-interval example below permit to visualize the numerical behavior of 
the algorithm. The folium example illustrates that, as expected, the non-convexity of 
K does not seem to penalize the moment approach. Finally, our experience reveals 
that the choice of alternative polynomial bases affects the quality of the approxima- 
tions. 

4.1. Single interval. Consider the elementary one-dimensional set K = [0, i] = 

{x G M : gi{x) — a;(i — a;) > 0} included in the unit interval B = [—1, 1]. We want 
to approximate vol(K) = i. Moments of the Lebesgue measure /12 on B are given 
byy2 = (2,0,2/3,0,2/5,0,2/7,...). 

Here is a simple Matlab script using GloptiPoly 3 instructions to input and solve 
the SDP relaxation Qd of the LP moment problem P with p = 1: 
>> d = 10; "/o degree 
>> mpol xO xl 

>> mO = meas(xO); ml = meas(xl); 
» gl = xl*(l/2-xl); 

» dm = (l+(0:d))'; y2 = ( (+1) . "dm-(-l) . ~dm) . /dm; 
>> yO = mom(mmon(xO,d) ) ; yl = mom(mmon(xl ,d) ) ; 

>> P = msdp(max(mass(ml) ) , gl>=0, y0==y2-yl) ; 7„ input moment problem 
>> msol(P); 7o solve SDP relaxation 
>> yl = double (mvec (ml) ) ; °/„ retrieve moment vector 

The volume estimate is then the first entry in vector yl. Note in particular the 
use of the moment constraint y0==y2-yl which ensures that moments yo of fiQ will 
be substituted by linear combinations of moments yi of /ii (decision variables) and 
moments y2 of /i2 (given) . 

Figure 14.11 displays three approximation sequences of vol (K) obtained by solving 
SDP relaxations (|3.1ip of increasing degrees d = 2, . . . , 50 of the infinite-dimensional 
LP moment problem P with three different parameters p: 
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Fig. 4.1. Three sequences of approximations o/vol[0, j] obtained by solving SDP relaxations 
of increasing degree. 



• the upper curve (in black) is a monotone non increasing sequence of upper 
bounds obtained by maximizing /d/ii, the mass of /ii, using the objective 
function max (mass (ml) ) in the above script; 

• the medium curve (in gray) is a sequence of approximations obtained by 
maximizing J pdiii with p := gi, using the objective function max(gl) in the 
above script; 

• the lower curve (in black) is a monotone non decreasing sequence of lower 
bounds on vol (K) obtained by computing upper bounds on the volume of 
B\K, using the objective function max (mass (ml) ) and the support constraint 
gl<=0 in the above script. The volume estimate is then 2-yl(l). 

We observe a much faster convergence when maximizing J gidfii instead of J dfii ; the 
upper and lower curves apparently exhibit slow convergence. 




Fig. 4.2. Positive polynomial approximation of degree 50 (solid) of the indicator function /,„ 
(dashed) on [—1,1]. 



To analyze these phenomena, we use solutions of the dual SDP problems, pro- 
vided automatically by the primal-dual interior-point method implemented in the 
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-Q.6 -0.4 -0.2 



Fig. 4.3. Positive polynomial approximation of degree 50 (solid) of the positive piecewise- 
polynomial function raa.y.{0,g\) on [—1,1], Polynomial gi is represented in dashed line. 




Fig. 4.4. Positive polynomial approximation of 
function 1 — /rg ii (dashed) on [—1, 1], 



50 (solid) of the complementary indicator 



SDP solver SeDuMi. On Figure W72\ we represent the degree-50 positive polynomial 
approximation h of the indicator function /k on B, which minimizes J„ hdx while 
satisfying /i — 1 > on K and ft, > on B \ K (yielding the volume estimate of the 
upper curve in Figure 14. ip . On Figure 14.31 we represent the degree-50 polynomial 
approximation h of the piecewise-polynomial function max(0,.gi), which minimizes 
Jg hdx while satisfying h — gi > on K and ft > on B \ K (yielding the volume es- 
timate of the medium curve in Figure WA} . On Figure H^ we represent the degree-50 
polynomial approximation ft of the complementary indicator function 1 — /k, which 
minimizes /g hdx while satisfying ft — 1 >OonB\K and ft > on K (yielding the 
volume estimate of the lower curve in Figure |4?T|) . We observe the characteristic oscil- 
lation phenomena near the boundary, typical of polynomial approximation problems 
|44j . The continuous function max(0, 51) is easier to approximate than discontinuous 
indicator functions, and this partly explains the better convergence of the medium 
approximation on Figure [¥?T] 

On Figures lT^ and H^ one observes relatively large oscillations near the boundary 
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points X £ { — 1,0, i,!} which significantly corrupt the quahty of the volume approx- 
imation. To some extent, these oscillations can be reduced by using a Chebyshev 
polynomial basis instead of the standard power basis. 




Fig. 4.5. Upper and lower bounds on vol [0. 
Chebyshev basis (black) and power basis (gray). 



obtained by solving SDP relaxations in the 



Figure 14.51 displays upper and lower bounds on the volume, computed up to 
degree 100, with the power basis (in gray) and with the Chebyshev basis (in black). 
Note that in order to input and solve SDP problems in the Chebyshev basis, we 
used our own implementation and the chebf un package since GloptiPoly 3 supports 
only the power basis. In Figure 14.51 we see that above degree 20 the quality of the 
bounds obtained with the power basis deteriorates, which suggests that the SDP solver 
encounters some numerical problems rather than convergence becoming slower (which 
is confirmed when changing to Chebyshev basis; see below). It seems that the SDP 
solver is not able to improve the bounds, most likely due to the symmetric Hankel 
structure of the moment matrices in the power basis: indeed, it is known that the 
conditioning (ratio of extreme singular values) of positive definite Hankel matrices is 
an exponential function of the matrix size jl9| . When the smallest singular values 
reach machine precision, the SDP solver is not able to optimize the objective function 
any further. 

In Figures 14.61 and 14.71 one observes that the degree-100 polynomial approxima- 
tion h{x) of the indicator function and its complement are tighter in the Cheby- 
shev basis (black) than in the power basis (gray). Firstly, we observe that the 
degree-100 approximations in the power basis do not significantly differ from the 
degree-50 approximations in the same basis, represented in Figures l4?2l and [44l This 
is consistent with the very fiat behavior of the right half of the upper and lower 
curves (in gray) in Figure H751 Secondly, some coefficients of h{x) in the power ba- 
sis have large magnitude h{x) = 1.0019 + 3.6161a; - 29.948.t2 -\ h SSUix'^^ + 

54985a;5° + . . . _ I018.4x^^ + 266692;^°° with the Euclidean norm of the coefficient 
vector greater than 10^. In contrast, the polynomial h(x) obtained in the Chebyshev 
basis h{x) = 0.1862io(a;) + 0.093432ti(x) - 0.30222i2(a;) + • • • + 0.0055367i49(a;) - 

0.020488^50(2:) + 0.0012267i99(a;) + 0.0011190tioo(a;) has a coefficient vector of 

Euclidean norm around 0.57627, where ifc(x) denotes the fc-th Chebyshev polynomial. 
Thirdly, oscillations around points x — and a; = 1/2 did not disappear with the 
Chebyshev basis, but the peaks are much thinner than with the power basis. Finally, 
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Fig. 4.6. Positive polynomial approximation of degree 100 of the indicator function I,^ i-, in 
the Chebyshev basis (black) and power basis (gray). 
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Fig. 4.7. Positive polynomial approximation of degree 100 of the indicator function 1 — .^rn 
in the Chebyshev basis (black) and power basis (gray). 



the oscillations near the interval ends x = —1 and x = 1 are almost suppressed, a 
well-known property of Chebyshev polynomials which have a denser root distribution 
near the interval ends. 

From these simple observations, we conjecture that a polynomial basis with a 
dense root distribution near the boundary of the semi-algebraic sets K and B should 
ensure a better convergence of the hierarchy of volume estimates. 

Finally, Figure 14.81 displays the CPU time required to solve the SDP problems 
(with SeDuMi, in the power basis in gray and in the Chebyshev basis in black) as a 
function of the degree, showing a polynomial dependence slightly slower than cubic 
in the power basis (due to the sparsity of moment matrices) and slightly faster than 
cubic in the Chebyshev basis. For example, solving the SDP problem of degree 100 
takes about 2.5 seconds of CPU time on our standard desktop computer. 

4.2. Bean. Consider K = {x G R^ : gi{x) = xi{xl+xl)-{x\ + x\xl'^xf) >{)} 
displayed in Figure HT91 which is a surface delimited by an algebraic curve gi{x) = 
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Fig. 4.8. CPU time required to solve the SDP relaxations (Chebyshev basis in black, power 
basis in gray) as a function of the degree. 




Fig. 4.9. Bean surface. 



of genus zero, hence rationally parametrizable. From the parametrization xi(t) = 
(1 + i^)/(l + t^ + t"^), X2{t) — txi{t), t G M, obtained with the algcurves package of 
Maple, we can calculate 



vol (K) 



- r Hz 



= /j^ dxidx2 = /r Xi{t)dx2{t) = /, 

- ZvIzL ^ 1.0581 



-t)(l+t)(l+t^)(l+3t'^+t"^) ^^ 



36 



with the help of the int integration routine of Maple. Similarly, we can calculate 
symbolically the first moments of the Lebesgue measure /ii on K, namely jjiqq = 
vol(K), yiio = ivol(K), yi^i - 0, yi^o = ivol(K), y^ = 0, yio2 = 3i^vol(K) 
etc. Observe that K C B = [-1, 1]^. 

On Figure 14.101 we represent a degree-20 positive polynomial approximation h 
of the indicator function In on B obtained by solving an SDP problem with 231 
unknown moments. We observe the typical oscillations near the boundary regions, 
but we can recognize the shape of Figure 14.91 

In Table H?T] we give relative errors in percentage observed when solving successive 
SDP relaxations (in the power basis) of the LP moment problems of maximizing 
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Fig. 4.10. Positive polynomial approximation of degree 20 of the indicator function of the bean 
surface. 



degree 


2 


4 


6 


8 


10 


12 


14 




error 


78% 


63% 


13% 


0.83% 


9.1% 


0.80% 


3.31% 




degree 


16 


18 


20 


22 


24 


26 


28 


30 


error 


3.8% 


3.3% 


2.6% 


5.6% 


4.1% 


4.1% 


3.9% 


3.7% 



Table 4.1 
Relative error when approximating the volume of the bean surface, as a function of the degree 
of the SDP relaxation. 



J gidfti. Note that the error sequence is not monotonically decreasing since we do 
not maximize J dfii and a good approximation can be obtained with few moments. 
Above degree 16, the approximation stagnates around 4%, Most hkely this is due to 
the use of the power basis, as already observed in the previous univariate examples. 
For example, at degree 20, one obtains the 6 first moment approximation 



y^fo - 1-10, y2To 



0.589, j/20? = 0.00, 2/220 = 0.390, j/211 = 0.00, j/202 



20 _ nnn „,_20 _ q -^22 



to be compared with the exact numerical values 

2/200 = 1-06, 2/210 = 0.579, 2/201 = 0.00, 2/220 = 0.386, 2/211 - 0.00, 2/202 = 0.119. 

Increasing the degree does not provide a better approximation. It is expected that a 
change of basis (e.g. multivariate Chebyshev or trigonometric) can be useful in this 
context. 

4.3. Folium. Consider K = {x e M^ : gi{x) = -{xf + x^f + \x\x\ > 0} 
displayed in Figure 14.111 which is a surface delimited by an algebraic curve of polar 
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Fig. 4.11. Folium surface. 



equation p — sm(26). The surface is contained in the unit disk B, on which the 
Lebesgue measure has moments 



y2a = 



(1 + (-i)"0(i + (-i)"-)r(i(i + ai))r(i(i + a,)) 

r(i(4 + ai+a2)) 



VaeN^ 



where F denotes the gamma function. The area is vol (K) — -^ Jq^ sin^{29)d9 — ^n 
and so, vol (K \ B) = tt - vol (K) = iyr. 

In Table l4?2l we give relative errors in percentage observed when solving successive 
SDP relaxations (in the power basis) of the LP moment problems of maximizing 
J gidfxi. We observe that nonconvexity of K does not play any special role. The 
quality of estimates does not really improve for degrees greater than 20. Here too, an 
alternative polynomial basis with dense root distribution near the boundaries of K 
and B would certainly help. 



degree 


4 


6 


8 


10 


12 


14 


16 


error 


87% 


19% 


14% 


9.4% 


4.3% 


4.5% 


5.9% 


degree 


18 


20 


22 


24 


26 


28 


30 


error 


1.2% 


5.3% 


5.9% 


7.2% 


8.7% 


9.0% 


8.8% 



Table 4.2 
Relative error when approximating the volume of the folium surface, as a function of the degree 
of the SDP relaxation. 



Figure 14.121 displays a degree-20 positive polynomial approximation h of the in- 
dicator function /k on B obtained by solving an SDP problem with 231 unknown 
moments. For visualization purposes, max(5/4, h) rather than h is displayed. Again 
typical oscillations occur near the boundary regions, but we can recognize the shape 
of Figure 14.111 
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■1 -' ^-^] 

Fig. 4.12. Positive polynomial approximation of degree 20 of the indicator function of the 
folium surface. 



5. Concluding Remarks. The methodology presented in this paper is general 
enough and applies to compact basic semi-algebraic sets which are neither necessarily 
convex nor connected. Its efficiency is related to the degree needed to obtain a good 
polynomial approximation of the indicator function of K (on a simple set that contains 
K) and from this viewpoint, convexity of K does not help much. On the other hand, 
the method is limited by the size of problems that SDP solvers presently available can 
handle. Moreover, the impact of the choice of the polynomial basis (e.g., Chebyshev or 
trigonometric) on the quality of the solution of the SDP relaxations deserves further 
investigation for a better understanding. Therefore, in view of the present status of 
SDP solvers and since in general high accuracy will require high degree, the method 
can provide good approximations for problems of small dimension (typically n = 2 
or n — 3). However, if one is satisfied with cruder bounds then one may consider 
problems in higher dimensions. 
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